home *** CD-ROM | disk | FTP | other *** search
/ Linux Cubed Series 8: LINUX Games / Linux Cubed Series 8 - LINUX Games.iso / games / video / fly8111-.000 / fly8111- / fly8 / oyplane.c < prev    next >
C/C++ Source or Header  |  1979-12-31  |  16KB  |  652 lines
  1. /* --------------------------------- oyplane.c ------------------------------ */
  2.  
  3. /* This is part of the flight simulator 'fly8'.
  4.  * Author: Eyal Lebedinsky (eyal@ise.canberra.edu.au).
  5. */
  6.  
  7. /* Dynamics of the second eXpermental plane.
  8.  * This model will eventualy do a proper aerodynamic simulation. It will
  9.  * go as far as 6DOF (simplified) and will properly implement a linearized
  10.  * serodynamics model. This program is being modified regularly as I learn
  11.  * more about the subject and gather more data. If I get a stable version
  12.  * then it will be separated into it's own module so as not to constantly
  13.  * have a half working program. But not yet.
  14.  *
  15.  * ToDo:
  16.  * 1 document units/scale used for all variables.
  17.  * 2 select aero formulas and break into classes:
  18.  *   - input parameters (from .prm)
  19.  *   - calculated fixed parameters
  20.  *   - slow changing parameters
  21.  *   - dynamically calculated variables.
  22.  * 3 use proper trans- and super-sonic calcs for Cd (wave drag).
  23.  * 4 add wing sweep to calcs?
  24. */
  25.  
  26. #include "plane.h"
  27.  
  28.  
  29. #define DOPITCH        (st.debug & DF_GPW)
  30.  
  31. /* 'CLf' is to avoid overflow (fractions are limited to the range [-2...+2)).
  32. */
  33. #define    CLf    4
  34.  
  35. LOCAL_FUNC void NEAR
  36. GetFoil (OBJECT *p, ANGLE flaps, ANGLE leFlaps, ANGLE *a0, int *maxCl,
  37.     int *minCl)
  38. {
  39.     int    t;
  40.  
  41.     *a0 = EP->Cl0 + fmul (EP->FEff, flaps);
  42.     t =  muldiv (EP->FEffCl, flaps, CLf);
  43.     t += muldiv (EP->lefEffCl, leFlaps, CLf);
  44.     *maxCl = EP->maxCl/CLf*10 + t;
  45.     *minCl = EP->minCl/CLf*10 + t;
  46. }
  47.     
  48. LOCAL_FUNC int NEAR
  49. GetCl (OBJECT *p, ANGLE aoa, ANGLE a0, int a, int maxCl, int minCl,
  50.     int *Cl, int *Fstall)
  51. {
  52.     int    aoaPstall, aoaNstall;
  53.     int    fstall, cl, stall;
  54.     ANGLE    t;
  55.  
  56.     aoa -= a0;                /* get absolute aoa */
  57.     a *= CLf;
  58.     aoaPstall = fmul (maxCl, a);
  59.     aoaNstall = fmul (minCl, a);
  60.     if (aoa > aoaPstall) {
  61.         cl = maxCl;
  62.         if (!(EX->flags & PF_NOSTALL)) {
  63.             t = aoa - aoaPstall;
  64.             if (t >= aoaPstall/2)
  65.                 fstall = 0;
  66.             else
  67.                 fstall = fdiv (t, aoaPstall)*2;
  68.         } else
  69.             fstall = FONE;
  70.         cl = fmul (cl, fstall);
  71.         stall = 1;
  72.     } else if (aoa < aoaNstall) {
  73.         cl = minCl;
  74.         if (!(EX->flags & PF_NOSTALL)) {
  75.             t = aoa - aoaNstall;
  76.             if (t <= aoaNstall/2)
  77.                 fstall = 0;
  78.             else
  79.                 fstall = fdiv (t, aoaNstall)*2;
  80.         } else
  81.             fstall = FONE;
  82.         cl = fmul (cl, fstall);
  83.         stall = 1;
  84.     } else {
  85.         fstall = FONE;
  86.         cl = fdiv (aoa, a);
  87.         stall = 0;
  88.     }
  89.     *Cl = cl;
  90.     if (Fstall)
  91.         *Fstall = fstall;
  92.     return (stall);
  93. }
  94.  
  95. /* We use SI units (Newton, Kilogram, meter, second, radian) but read in some
  96.  * data in other units (feet, pound, pound-force, degrees etc.).
  97.  *
  98.  * Since we use integers, some scaling factors are introduced, sometimes
  99.  * varying along the calculations. Angles are stored in a special format
  100.  * were pi/2 radians (90 degrees) is 1.0, and the representation wraps around
  101.  * at 2.0, the full range being [-2.0...2.0).
  102.  *
  103.  * Definitions:
  104.  *
  105.  * b    wing span (m)
  106.  * S    wing area (m^2)
  107.  * ClRate0 slope of Cl for the wing foil (1/r)
  108.  *    we assume it to be 2*pi 
  109.  * e    Oswald's (or span) efficiency factor
  110.  * aoa0    angle of attack where Cl is zero (r)
  111.  * feff    degrees aoa effective for 1 degree flaps
  112.  *
  113.  * AR    wing aspect ratio
  114.  * k_AR    1/(pi*e*AR)
  115.  * ClRate  slope of Cl for wing (1/r)
  116.  * e1    drag span-efficiency-factor. We use 'e' for it.
  117.  *
  118.  * V    airspeed (m/s)
  119.  * rho    air density (k/m^3)
  120.  * aoa    angle of attack, geometric (r)
  121.  * aoaf    angle of attack effect due to flaps (r)
  122.  * aoai    induced angle of attack (r)
  123.  *    = Cl / (pi*e1*AR), but folded into ClRate
  124.  * aeff    effective angle of attack (r)
  125.  *  ... all of these angles in degrees!
  126.  * q    dynamic pressure (k/m^2)
  127.  * Cl    lift coefficient
  128.  * Cd    drag coefficient
  129.  * Cdi    induced Cd
  130.  * Cdp    profile Cd
  131.  *
  132.  * Formulas used:
  133.  *
  134.  * AR = b^2 / S
  135.  * aoai = Cl / (pi*e1*AR)
  136.  * ClRate = ClRate0 / (1+ClRate0/(pi*e1*AR))
  137.  *  or:
  138.  * 1/ClRate = 1/ClRate0 + 1/(pi*e1*AR)
  139.  *  since we assume ClRate0 = 2*pi we now have
  140.  * 1/ClRate = 2*pi*e1*AR/(e1*AR+2)
  141.  *
  142.  * aoaf = flaps * feff
  143.  * aeff = aoa + aoaf
  144.  * Cl = (aeff - aoa0) * ClRate
  145.  *
  146.  *
  147.  *
  148.  *
  149.  * Cdi = Cl^2 /(pi*e*AR)
  150.  * Cdp = Cdp_characteristic + Cdp_flaps + Cdp_gear + Cdp_bombs + Cdp_speedbrakes
  151.  * Cd = Cdi + Cdp
  152.  * q = rho * V^2 / 2
  153.  * lift is perpendicular to airflow:
  154.  *     lift = Cl * q * S
  155.  * drag is parallel to airflow:
  156.  *     drag = Cd * q * S
  157.  *
  158.  * We need these factors:
  159.  * Ixx    pitching moment of inertia
  160.  * Iyy    rolling moment of inertia
  161.  * Izz    yawing moment of inertia
  162.  * we calculate the moments as:
  163.  * Mx = lift*Lc + drag*Dc + thrust*Tc + (stablilizer+elevators)*Sc
  164.  * Lc, Dc, Tc, Sc: effective moments levers.
  165.  * My,Mz = f(many things...)
  166. */
  167.  
  168. extern void FAR
  169. dynamics_yplane (OBJECT *p, int action)
  170. {
  171.     POINTER    *ptr;
  172.     int    onground;
  173.     int    t, weight, thrust, wfact;
  174.     int    rrho, rho, sos, sa, ca, sb, cb;
  175.     int    stall;
  176.     long    tt;
  177.  
  178.     if (action)
  179.         return;
  180.  
  181.     if (dynamics_input (p))
  182.         return;
  183.  
  184.     ptr = p->pointer;
  185.  
  186.     onground = EX->flags & PF_ONGROUND;
  187.  
  188.     if (onground ? check_takeoff (p) : check_land (p)) {
  189.         p->flags |= F_HIT;
  190.         return;
  191.     }
  192.  
  193. /* Automatically refuel when stationery on ground.
  194. */
  195.     if (onground && 0 == p->speed && 0 == EX->power)
  196.         supply (p, 0);
  197.  
  198.     airdata (p->R[Z], 0, &rrho, &rho, &sos);
  199. CFshow (p, "rrho", rrho);
  200. CFshow (p, "rho", rho);
  201. CCshow (p, 0, "sos", (long)sos/VONE);
  202.  
  203. /* Get engine thrust.
  204. */
  205.     f16engine (p, sos);
  206.  
  207. /* We select a factor to use in the integer arithmetic to avoid overflow
  208.  * and maintain resolution.
  209. */
  210.     wfact = (int)(EP->weight/1000)+1;
  211. CCshow (p, 0, "wfact", (long)wfact);
  212.  
  213. /* Thrust is stored as lbf/10.
  214. */
  215.     t = muldiv (4*10, EX->thrust, wfact);
  216.     thrust = fmul (t, FCON(4.4482/4));            /* lbf->N */
  217. CCshow (p, 0, "athrust", (long)thrust*wfact);
  218.  
  219.     tt = EP->weight + EX->fuel/100;                /* lb */
  220.     if ((t = EX->stores[WE_M61-1]) > 0)
  221.         tt += t * st.bodies[O_M61]->shape->weight/1000;
  222.     if ((t = EX->stores[WE_MK82-1]) > 0)
  223.         tt += t * st.bodies[O_MK82]->shape->weight/1000;
  224.     weight = fmul ((int)(tt/wfact), FCON (0.4536));        /* Kg */
  225. CCshow (p, 0, "weight", (long)weight*wfact);
  226.  
  227.     thrust = muldiv (thrust, VONE, weight);
  228. CVshow (p, "thrust", thrust);
  229.  
  230. /* Flying, although may still have ground contact. Must do full aerodynamics.
  231. */
  232.     {
  233.     typedef int    FRACTION;
  234.  
  235.     FRACTION    AR;        /* wing aspect ratio (1/AR) */
  236.     FRACTION    pi_e_AR;    /* 1/(e*pi*AR) */
  237.     FRACTION    a;        /* Cl rate */
  238.     ANGLE        aFlaps;        /* flaps */
  239.     ANGLE        leFlaps;    /* leading edge flaps */
  240.     ANGLE        alpha;        /* angle of attack (aoa) */
  241.     ANGLE        beta;        /* angle of sideslip */
  242.     FRACTION    Cl;
  243.     int        Fstall;
  244.     FRACTION    Cdi;
  245.     FRACTION    Cd;
  246.     int        S;
  247.     int        b;
  248.     int        c;
  249.     int        V;
  250.     int        qS_V;
  251.     int        qS;
  252.     int        qS_V_b2;
  253.     int        qS_V_c;
  254.     int        lift;
  255.     int        drag;
  256.     int        sf;        /* side force */
  257.     int        b2fact;
  258.     int        qSfact;
  259.     int        m;
  260.     FRACTION    Crudder, Cailerons, Celevators;
  261.     FRACTION    Ca, Cb, Cdx, Cdy, Cdz;
  262.     VECT        F, M;
  263.     MAT        I;
  264.  
  265. /* First get some basic data.
  266. */
  267.     AR = fdiv (EP->wing_area, muldiv (EP->wing_span, EP->wing_span, VONE));
  268.                         /* AR */
  269.     t = fmul (AR, FCON(1.0/C_PI));        /* 1/(pi*AR) */
  270.     pi_e_AR = muldiv (t, 100, EP->efficiency_factor);
  271.                         /* 1/(e*pi*AR) */
  272.  
  273.     a = FCON(1/(2*C_PI));            /* note we keep 1/a */
  274.     a += pi_e_AR;
  275.     a = RAD2ANG (a);
  276.  
  277. /* Find our alpha, beta and friends.
  278. */
  279.     if (onground && iabs(EX->v[Y]) < 2*VONE)
  280.         alpha = beta = 0;
  281.     else {
  282.         alpha = -ATAN (EX->v[Z], EX->v[Y]);
  283.         beta  =  ASIN (fdiv (EX->v[X], p->speed));
  284.     }
  285.     alpha += EP->Aoffset;            /* wing rigging offset */
  286.     EX->aoa = alpha;
  287.     sa = SIN (alpha);
  288.     ca = COS (alpha);
  289.     sb = SIN (beta);
  290.     cb = COS (beta);
  291. CAshow (p, "alpha", alpha);
  292. CAshow (p, "beta", beta);
  293.  
  294. /* Get flaps effect. We introduce automatic flaps in response to alpha.
  295. */
  296.     aFlaps = EX->flaps;
  297.     if (EX->flags & PF_AUTOFLAP) {
  298.         if (aFlaps) {
  299.             EX->leFlaps = aFlaps;
  300.         } else {
  301.             if (alpha > EP->AFamin) {
  302.                 aFlaps = fmul (alpha-EP->AFamin,
  303.                                 EP->AFrate);
  304.                 if (aFlaps > EP->AFmax)
  305.                     aFlaps = EP->AFmax;
  306.             } 
  307.             if (alpha > EP->ALEFamin) {
  308.                 EX->leFlaps = fmul (alpha-EP->ALEFamin,
  309.                                 EP->ALEFrate);
  310.                 if (EX->leFlaps > 100)
  311.                     EX->leFlaps = 100;
  312.             } else
  313.                 EX->leFlaps = 0;
  314.         }
  315.     } else {
  316.         EX->leFlaps = 0;
  317.     }
  318.     aFlaps  = muldiv (EP->MaxFlaps, aFlaps, 100);
  319.     leFlaps = muldiv (EP->MaxLEFlaps, EX->leFlaps, 100);
  320. {
  321.     int    maxCl, minCl;
  322.     ANGLE    a0;
  323.  
  324.     GetFoil (p, aFlaps, leFlaps, &a0, &maxCl, &minCl);
  325.     stall = 0;
  326.     EX->flags &= ~PF_STALL;
  327.     if (GetCl (p, alpha, a0, a, maxCl, minCl, &Cl, &Fstall)) {
  328.         EX->flags |= PF_STALL;
  329.         if (!(EX->flags & PF_NOSTALL))
  330.             stall = 1;
  331.     }
  332. }
  333.     Cdi = fmul (pi_e_AR, Cl)*CLf;
  334.     Cdi = fmul (Cdi, Cl)*CLf;
  335.  
  336. CFshow (p, "Cl", Cl);
  337. CFshow (p, "Cdi", Cdi);
  338.  
  339. /* The groung effect formula is extracted from the graph in Smiths 'The
  340.  * Illustrated Guide To Aerodynamics' end of chapter 3.
  341. */
  342.     t = EP->wing_span;
  343.     if (p->R[Z] < (long)t) {        /* ground effect */
  344.         t = fdiv ((int)p->R[Z], t);    /* h/b */
  345. CFshow (p, "h/b", t);
  346.         if (t < FCON(0.1))
  347.             t = 5 * t;
  348.         else
  349.             t = FCON(1.06) - fdiv (FCON(0.07), t);
  350. CFshow (p, "gef", t);
  351.         Cdi = fmul (Cdi, t);
  352. CFshow (p, "Cdi", Cdi);
  353.     }
  354.  
  355.     Cd = EP->Cdp0;
  356. CFshow (p, "Cdp0", Cd);
  357.     Cd += Cdi;
  358. CFshow (p, "+Cdi", Cd);
  359.     if (EX->airbrake) {
  360.         Cd += muldiv (EP->Cds, EX->airbrake, 100);
  361. CFshow (p, "+Brks", Cd);
  362.     }
  363.     if (EX->equip & EQ_GEAR) {
  364.         Cd += EP->Cdg;
  365. CFshow (p, "+Gear", Cd);
  366.     }
  367.     if (EX->stores[WE_MK82-1] > 0) {
  368.         Cd += EX->stores[WE_MK82-1] * EP->CdMK82;
  369. CFshow (p, "+MK82", Cd);
  370.     }
  371.  
  372. /* We now calculate the new forces, based on the state of the plane.
  373.  * The given state is:
  374.  *  Cdx        Rotation around x (nose up, pitch)
  375.  *  Cdy        Rotation around y (right wing down, roll)
  376.  *  Cdz        Rotation around z (nose right, yaw)
  377.  *  Crudder    rudder left-turn angle
  378.  *  Cailerons    ailerons roll-right (cwise) angle
  379.  *  Celevators    elevators push-down angle
  380.  *  Ca        alpha: nose above wind angle
  381.  *  Cb        beta:  nose left-of-wind angle
  382.  * The forces are:
  383.  *  Fx        right
  384.  *  Fy        forward
  385.  *  Fz        up
  386.  * The moments are:
  387.  *  Mx        around x (nose up, pitch)
  388.  *  My        around y (right wing down, roll)
  389.  *  Mz        around z (nose right, yaw)
  390.  *
  391.  * From each force we get acceleration [a = F/m] and then velocity [v += a*t]
  392.  * and position [x += (v + a*t/2)*t].
  393.  *
  394.  * From each moment we get angular acceleration [aa = M/I], then we get
  395.  * rotation rate [da += aa*t] and angle [a += (da + aa*t/2)*t].
  396. */
  397.     V = p->speed/VONE;
  398. CCshow (p, 0, "V", (long)V);
  399.     S = EP->wing_area;
  400. CCshow (p, 0, "S", (long)S/VONE);
  401.     b = EP->wing_span;
  402. CCshow (p, 2, "b", b*100L/VONE);
  403.     c = EP->wing_cord;
  404. CCshow (p, 2, "c", c*100L/VONE);
  405.  
  406. /* Scale factors: These were selected so avoid integer overflow while
  407.  * maintaining good resolution, so they may look weird.
  408.  *
  409.  * item:    is larger than real life
  410.  *        by a factor of:
  411.  *
  412.  * b        *VONE
  413.  * c        *VONE
  414.  * S        *VONE
  415.  * qS_V        *VONE/wfact
  416.  * qS        *VONE/W
  417.  * qS_V_b2    /VONE/W            [frac]
  418.  * qS_V_c    /VONE/W            [frac]
  419.  * Cdx        *2/PI/VONE        [frac]
  420.  * Cdy        *2/PI/VONE        [frac]
  421.  * Cdz        *2/PI/VONE        [frac]
  422.  * lift        *VONE/W
  423.  * drag        *VONE/W
  424.  * side force    *VONE/W
  425.  * F[]        *VONE/W
  426.  * M[]        *2/PI/VONE/VONE/W    [frac]
  427.  * I[]        /W/10            [frac]
  428. */
  429.     qS_V    = muldiv (fmul (rho, S), V, 2*wfact);
  430.     qS      = muldiv (qS_V, V, weight);
  431.     qS_V_b2 = muldiv (qS_V, (FONE/VONE/VONE/VONE)*b/2, weight);
  432.     qS_V_c  = muldiv (qS_V, (FONE/VONE/VONE/VONE)*c, weight);
  433. CCshow (p, 0, "qS", (long)qS);
  434. CCshow (p, 0, "qS_V_b2", (long)qS_V_b2);
  435. CCshow (p, 0, "qS_V_c", (long)qS_V_c);
  436.  
  437. /* Keep the load within the designated range by limiting the elevators.
  438. */
  439.     Celevators = EX->elevators + EX->tElevators;
  440.     Celevators = fmul (Celevators, Fstall);
  441.     if (EX->flags & PF_AUTOELEVATOR) {
  442.         t = EP->AErate/VONE;
  443.         if (V > t)
  444.             Celevators = muldiv (Celevators, t, V);
  445.     }
  446.     Celevators = -muldiv (EP->MaxElevators, Celevators, RAD2ANG(100));
  447.  
  448.     Cailerons = -muldiv (EP->MaxAilerons, EX->ailerons, RAD2ANG(100));
  449.     Cailerons = fmul (Cailerons, Fstall);
  450.  
  451.     Crudder = muldiv (EP->MaxRudder, EX->rudder+EX->tRudder, RAD2ANG(100));
  452.     Crudder = fmul (Crudder, Fstall);
  453.  
  454.     Cdx =  p->da[X];
  455.     Cdy =  p->da[Y];
  456.     Cdz = -p->da[Z];
  457.  
  458.     if (alpha >= FCON(4/C_PI))
  459.         Ca = FCON(1.99);
  460.     else if (alpha <= -FCON(4/C_PI))
  461.         Ca = -FCON(1.99);
  462.     else
  463.         Ca = ANG2RAD(alpha);
  464.  
  465.     Cb = ANG2RAD(beta);
  466. CFshow (p, "Crdr", Crudder);
  467. CFshow (p, "Cail", Cailerons);
  468. CFshow (p, "Celev", Celevators);
  469. CFshow (p, "Cdx", Cdx);
  470. CFshow (p, "Cdy", Cdy);
  471. CFshow (p, "Cdz", Cdz);
  472. CFshow (p, "Ca", Ca);
  473. CFshow (p, "Cb", Cb);
  474.  
  475.     M[X] = M[Y] = M[Z] = 0;
  476.  
  477. /* Lift, drag and side-force.
  478. */
  479.     lift =  fmul (Cl, qS * CLf);
  480.     drag = -fmul (Cd, qS);
  481.     sf   =    fmul (qS, fmul (EP->Cydr, Crudder));
  482.     if (Cb < FCON(0.1))
  483.         sf += fmul (qS, fmul (EP->Cybeta, Cb*10));
  484.     else {
  485.         t = fmul (EP->Cybeta, Cb);
  486.         if (t < FCON(0.2))
  487.             sf += fmul (qS, t*10);
  488.         else
  489.             sf += fmul (qS, t)*10;
  490.     }
  491. CVshow (p, "lift", lift);
  492. CVshow (p, "drag", drag);
  493. CVshow (p, "sf", sf);
  494.  
  495. /* Aerodynamic forces about the c.g..
  496. */
  497.     F[X] = fmul (sb, drag) + sf;
  498.     t    = fmul (cb, drag);
  499.     F[Y] = fmul (sa, lift) + fmul (ca, t);
  500.     F[Z] = fmul (ca, lift) - fmul (sa, t);
  501.  
  502. /* Engine forces about the c.g..
  503.  *
  504.  * We break the thrust into linear and angular components.
  505. */
  506.     F[Y] += fmul (thrust, COS(EP->Ea));
  507.     F[Z] += fmul (thrust, SIN(EP->Ea));
  508.  
  509.     t = -fmul (thrust, SIN(EP->Ea + EP->Eb));    /* angular */
  510.     M[X] += fmul (t, EP->Er);
  511.  
  512. CVshow (p, "Tx", F[X]);
  513. CVshow (p, "Ty", F[Y]);
  514. CVshow (p, "Tz", F[Z]);
  515.  
  516. /* Add force of gravity.
  517. */
  518.     F[X] -= fmul (GACC, p->T[X][Z]);
  519.     F[Y] -= fmul (GACC, p->T[Y][Z]);
  520.     F[Z] -= fmul (GACC, p->T[Z][Z]);
  521.  
  522. /* Attempt to improve resolution.
  523. */
  524. #define MAXB2    FCON(0.5)
  525.     if (qS_V_b2 < MAXB2/(VONE*VONE))    /* too small */
  526.         b2fact = VONE*VONE;
  527.     else if (F(b2fact = MAXB2/qS_V_b2))    /* just right */
  528.         b2fact = 1;             /* too large */
  529. CCshow (p, 0, "b2fact", b2fact);
  530.     qS_V_b2 *= b2fact;
  531.     qS_V_c  *= b2fact;
  532.     b2fact    *= VONE;
  533. #undef MAXB2
  534.  
  535. #define MAXQS    FCON(C_PI/2*VONE*VONE*VONE/FONE)
  536.     if (qS < MAXQS/(VONE*VONE))        /* too small */
  537.         qSfact = VONE*VONE;
  538.     else if (F(qSfact = MAXQS/qS))        /* just right */
  539.         qSfact = 1;             /* too large */
  540. CCshow (p, 0, "qSfact", qSfact);
  541.     qS      *= qSfact;
  542.     qSfact    *= VONE;
  543.  
  544.     qS = fdiv (qS, MAXQS);
  545. #undef MAXQS
  546.  
  547.  
  548. /* Pitching moment.
  549. */
  550.  
  551. if (DOPITCH) {
  552.     m =    + fmul (qS, EP->Cm0)
  553.         + fmul (fmul (qS, EP->Cmalpha), Ca);
  554.     M[X] += muldiv (m, c, qSfact);
  555. CFshow (p, "Cmx", m);
  556. }
  557.     m =    + fmul (fmul (qS, EP->Cmde), Celevators);
  558.     M[X] += muldiv (m, c, qSfact);
  559.  
  560.     m =    + fmul (fmul (qS_V_c, EP->Cmq), Cdx);
  561.     M[X] += muldiv (m, 10*c, b2fact);
  562.  
  563.     
  564. /* Rolling moment.
  565. */
  566.     m =    + fmul (fmul (qS, EP->Clda), Cailerons)
  567.         + fmul (fmul (qS, EP->Cldr), Crudder)
  568.         + fmul (fmul (qS, EP->Clbeta), Cb);    /* hihedral */
  569.     M[Y] += muldiv (m, b, qSfact);
  570.  
  571. #define Clr    (Cl*(CLf/4))
  572.     m =    + fmul (fmul (qS_V_b2, EP->Clp), Cdy)    /* damping */
  573.         + fmul (fmul (qS_V_b2, Clr), Cdz);    /* yaw induced */
  574.     M[Y] += muldiv (m, b, b2fact);
  575. #undef Clr
  576.  
  577.  
  578. /* Yawing moment.
  579. */
  580.     m =    + fmul (fmul (qS, EP->Cndr), Crudder)
  581.         + fmul (fmul (qS, EP->Cnda), Cailerons)
  582.         + fmul (fmul (qS, EP->Cnbeta), Cb);    /* weathercock */
  583.     M[Z] -= muldiv (m, b, qSfact);
  584.  
  585.     m =    + fmul (fmul (qS_V_b2, EP->Cnr), Cdz)    /* damping */
  586.         + fmul (fmul (qS_V_b2, EP->Cnp), Cdy);
  587.     M[Z] -= muldiv (m, b, b2fact);
  588.  
  589. CFshow (p, "Mx", M[X]);
  590. CFshow (p, "My", M[Y]);
  591. CFshow (p, "Mz", M[Z]);
  592.  
  593.     I[X][X] = EP->Iyy;
  594.     I[Y][Y] = EP->Ixx;
  595.     I[Z][Z] = EP->Izz;
  596.     I[Y][Z] = EP->Izx;
  597.  
  598.     LandGear (p, F, M);
  599.  
  600.     SixDOF (p, F, M, I);
  601.  
  602.     Mobj (p);
  603.  
  604. CCshow (p, 3, "G", EX->Gforce*1000L/GACC);
  605. tt = p->R[Z] + V*VONE*(long)V/(2L*GACC/VONE);
  606. CCshow (p, 0, "He", (long)(tt*3.28/VONE));
  607. tt = p->V[Z] + VONE*V*(long)ihypot3d(EX->a)/GACC;
  608. CCshow (p, 0, "Ps", (long)(tt*3.28/VONE));
  609. tt = tt*60L;
  610. CCshow (p, 0, "RC", (long)(tt*3.28/VONE));
  611. }
  612.  
  613.     p->speed = ihypot3d (EX->v);
  614.     VMmul (p->V, EX->v, p->T);
  615.  
  616. #define MAX_SPEED    (1000*VONE)
  617.     if (p->speed > MAX_SPEED) {            /* temp */
  618.         t = fdiv (MAX_SPEED, p->speed);
  619.         EX->v[X] = fmul (t, EX->v[X]);
  620.         EX->v[Y] = fmul (t, EX->v[Y]);
  621.         EX->v[Z] = fmul (t, EX->v[Z]);
  622.         p->V[X]  = fmul (t, p->V[X]);
  623.         p->V[Y]  = fmul (t, p->V[Y]);
  624.         p->V[Z]  = fmul (t, p->V[Z]);
  625.         p->speed = ihypot3d (EX->v);
  626.     }
  627. #undef MAX_SPEED
  628.  
  629.     if (onground)
  630.         EX->flags &= ~PF_STALL;
  631.  
  632.     if (EX->Gforce > EP->max_lift || EX->Gforce < EP->min_lift)
  633.         EX->flags |= PF_GLIMIT;
  634.     else
  635.         EX->flags &= ~PF_GLIMIT;
  636.  
  637. /* Mach number.
  638. */
  639.     EX->mach = muldiv (p->speed, 1000, sos);
  640.  
  641. /* pull up warning time.
  642. */
  643.     t = muldiv (10000, iabs(p->speed), 300*VONE);
  644.     t = muldiv (t, iabs(p->a[X]), D90);
  645.     if (t < 2000)
  646.         t = 2000;
  647.     EX->misc[8] = t;
  648. }
  649.  
  650. #undef DOPITCH
  651. #undef CLf
  652.